General Model Approach A - Univariate
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(scoringRules)
library(ggplot2)
# Load the model data
# Here: NeuralProphet
model <- read.csv("../data/forecasts/peaks_22-24_model-neuralprophet.csv")
model$ds <- as.POSIXct(model$ds, tz = "UTC")
# Here: Arma
# model <- read.csv("./data/forecasts/peaks_22-24_model-arma.csv")
# model$ds <- as.POSIXct(model$ds, tz = "UTC")
# specify the length for the error learning phase in days
length <- 365
# Filter only for the years 2023 (Training) and 2024 (Testing)
model <- model |>
filter((year(ds) == 2023) | (year(ds) == 2024))
getCRPS_A <- function(d) {
# Parameter: d day to predict
print(d)
# --------- Error Learning Phase ---------
# Getting the test data: starting from day i get the next 365 days
# Achtung: Hier dürfen nur Training oder Testing data verwendet werden
df_train <- model[d:(364 + d), ]
# Standardize the residuals
mu_resid <- mean(df_train$residuals)
sigma_resid <- sqrt(var(df_train$residuals))
df_train$residuals_std <- (df_train$residuals - mu_resid) / sigma_resid
# hist(df_train$residuals_std, main = "Histogram of Standardized Residuals", xlab = "Standardized Residuals")
# Generate the Empirical Distribution Function for the resids
ecdf <- ecdf(df_train$residuals_std)
# Define the inverse ECDF
# Input p is the probability and the return is the quantile
inverse_ecdf <- function(p) {
quantile(ecdf, p, names = FALSE)
}
print("Error Learning DONE...")
# ------- Prediction phase -------
m <- 90
# Define the quantiles
quantiles <- seq(1 / (m + 1), m / (m + 1), 1 / (m + 1))
peaks_dis_A <- data.frame(quantiles = quantiles, values = rep(0, length(quantiles)))
# Get the next 24 hours after the testing period
# Hier nur Testing data vom Model!!!
df_test <- model[(d + length), ]
for (i in quantiles) {
# Destandardize the error for the qunatile i
quantile_resid <- (inverse_ecdf(i) * sigma_resid) + mu_resid
peaks_dis_A[peaks_dis_A$quantiles == i, "values"] <- df_test$yhat[1] + quantile_resid
}
print("Prediction DONE...")
# Histogram
hist(peaks_dis_A[, "values"], main = "Histogram of Peaks Distribution A", xlab = "Peaks", breaks = "Sturges")
# Extracting the peak from the test day
peak <- df_test[, "y"]
# CRPS-Score
crps_sample(peak, peaks_dis_A[, "values"])
}
# specify the length for rolling iterations in days
len_test <- 28
crps_scores <- list()
for (d in seq(1, len_test)) {
crps_score <- getCRPS_A(d)
# Append the CRPS score to the list
crps_scores[[d]] <- crps_score
}
## [1] 1
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 2
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 3
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 4
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 5
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 6
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 7
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 8
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 9
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 10
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 11
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 12
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 13
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 14
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 15
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 16
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 17
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 18
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 19
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 20
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 21
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 22
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 23
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 24
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 25
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 26
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 27
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 28
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
print("Mean CRPS for 28 days in 2024 for NeuralProphet")
## [1] "Mean CRPS for 28 days in 2024 for NeuralProphet"
print(mean(unlist(crps_scores)))
## [1] 2777.669
library(tidyverse)
library(scoringRules)
library(ggplot2)
# Load the model data
# Here: NeuralProphet
# model <- read.csv("./data/forecasts/peaks_22-24_model-neuralprophet.csv")
# model$ds <- as.POSIXct(model$ds, tz = "UTC")
# Here: Arma
model <- read.csv("../data/forecasts/peaks_22-24_model-arma.csv")
model$ds <- as.POSIXct(model$ds, tz = "UTC")
# specify the length for the error learning phase in days
length <- 365
# Filter only for the years 2023 (Training) and 2024 (Testing)
model <- model |>
filter((year(ds) == 2023) | (year(ds) == 2024))
getCRPS_A <- function(d) {
# Parameter: d day to predict
print(d)
# --------- Error Learning Phase ---------
# Getting the test data: starting from day i get the next 365 days
# Achtung: Hier dürfen nur Training oder Testing data verwendet werden
df_train <- model[d:(364 + d), ]
# Standardize the residuals
mu_resid <- mean(df_train$residuals)
sigma_resid <- sqrt(var(df_train$residuals))
df_train$residuals_std <- (df_train$residuals - mu_resid) / sigma_resid
# hist(df_train$residuals_std, main = "Histogram of Standardized Residuals", xlab = "Standardized Residuals")
# Generate the Empirical Distribution Function for the resids
ecdf <- ecdf(df_train$residuals_std)
# Define the inverse ECDF
# Input p is the probability and the return is the quantile
inverse_ecdf <- function(p) {
quantile(ecdf, p, names = FALSE)
}
print("Error Learning DONE...")
# ------- Prediction phase -------
m <- 90
# Define the quantiles
quantiles <- seq(1 / (m + 1), m / (m + 1), 1 / (m + 1))
peaks_dis_A <- data.frame(quantiles = quantiles, values = rep(0, length(quantiles)))
# Get the next 24 hours after the testing period
# Hier nur Testing data vom Model!!!
df_test <- model[(d + length), ]
for (i in quantiles) {
# Destandardize the error for the qunatile i
quantile_resid <- (inverse_ecdf(i) * sigma_resid) + mu_resid
peaks_dis_A[peaks_dis_A$quantiles == i, "values"] <- df_test$yhat[1] + quantile_resid
}
print("Prediction DONE...")
# Histogram
hist(peaks_dis_A[, "values"], main = "Histogram of Peaks Distribution A", xlab = "Peaks", breaks = "Sturges")
# Extracting the peak from the test day
peak <- df_test[, "y"]
# CRPS-Score
crps_sample(peak, peaks_dis_A[, "values"])
}
# specify the length for rolling iterations in days
len_test <- 28
crps_scores <- list()
for (d in seq(1, len_test)) {
crps_score <- getCRPS_A(d)
# Append the CRPS score to the list
crps_scores[[d]] <- crps_score
}
## [1] 1
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 2
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 3
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 4
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 5
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 6
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 7
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 8
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 9
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 10
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 11
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 12
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 13
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 14
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 15
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 16
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 17
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 18
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 19
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 20
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 21
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 22
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 23
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 24
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 25
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 26
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 27
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
## [1] 28
## [1] "Error Learning DONE..."
## [1] "Prediction DONE..."
print("Mean CRPS for 28 days in 2024 for ARMA")
## [1] "Mean CRPS for 28 days in 2024 for ARMA"
print(mean(unlist(crps_scores)))
## [1] 7875.091
General Model Approach B - Multivariate
library(tidyverse)
library(scoringRules)
library(ggplot2)
library(copula)
##
## Attaching package: 'copula'
## The following object is masked from 'package:lubridate':
##
## interval
# This is a first general implementation of version B
# Load the model data
# Here: NeuralProphet
model <- read.csv("../data/forecasts/load_22-24_model-neuralprophet_2024IsForecasted.csv")
model$ds <- as.POSIXct(model$ds, tz = "UTC")
# Here: Arma
# model <- read.csv("./data/forecasts/load_22-24_model-arma.csv")
# model$ds <- as.POSIXct(model$ds, tz = "UTC")
# specify the length for the error learning phase in days
length <- 365
# Filter only for the years 2023 (Training) and 2024 (Testing)
model <- model |>
filter((year(ds) == 2023) | (year(ds) == 2024))
getCRPS_B <- function(d) {
print(d)
# --------- Error Learning Phase ---------
# Getting the test data: starting from day i get the next 365 days
# Achtung: Hier dürfen nur Training oder Testing data verwendet werden
df_train <- model[((d - 1) * 24 + 1):(((364 + d) * 24) - 1), ]
mu_sigma <- data.frame(hour = seq(0, 23), mu = rep(0, 24), sigma = rep(0, 24))
for (i in seq(0, 23)) {
load_h <- df_train |>
filter(hour(ds) == i)
mean <- mean(load_h$residuals)
std <- sqrt(var(load_h$residuals))
mu_sigma[i + 1, "mu"] <- mean
mu_sigma[i + 1, "sigma"] <- std
# Extract the residuals from the data
resid <- data.frame(resid = load_h$residuals)
# Standardize the residuals
resid <- resid |>
mutate(resid_std = (resid - mu_sigma[i + 1, "mu"]) / mu_sigma[i + 1, "sigma"])
# Generate the ECDF
assign(paste0("ecdf_", i), ecdf(resid$resid_std))
}
# Generate the Inverse ECDF
# Input: p... Probabilty , i... hour
# Output: quantile
inverse_ecdf <- function(p, i) {
quantile(get(paste0("ecdf_", i)), p, names = FALSE)
}
print("Error Learning DONE...")
# ---- Dependence Learning Phase ----
# Define the length m of the learning phase
m <- 90
# Initialize the training sample of the past m errors for every hour to learn the copula
X <- matrix(nrow = m, ncol = 0)
for (i in seq(0, 23)) {
# get the hour data
load_h <- df_train |>
filter(hour(ds) == i)
# get the errors
epsilon <- as.vector(load_h$residuals)
# Select only the last m errors
epsilon <- epsilon[(length(epsilon) - m):length(epsilon)]
# standardize the errors
epsilon_std <- (epsilon - mu_sigma[i + 1, "mu"]) / mu_sigma[i + 1, "sigma"]
# Apply the ecdf
ecdf <- get(paste0("ecdf_", i))
epsilon <- ecdf(epsilon_std)
X <- cbind(X, epsilon)
}
# Create the Copula and the Rank Matrix
# 2) Non-Parametric Approach
# Consider X as the sample from the empirical copula
# Apply the rank function to get the rank matrix
R_emp <- apply(X, 2, rank)
print("Dependence Learning DONE...")
# ------- Prediction phase -------
# Generate the 24 univariate for every hour
# m also defines the quantile level i.. 1 - m with i/m+1 quantiles
univariate_forecast_t <- matrix(nrow = 24, ncol = m)
quantiles <- seq(1 / (m + 1), m / (m + 1), 1 / (m + 1))
# Get the next 24 hours after the testing period
# Hier nur Testing data vom Model!!!
df_test <- model[((d + length - 1) * 24):((d + length) * 24 - 1), ]
predict_point <- df_test[, "yhat"]
for (i in seq(0, 23)) {
# Get the prediction for the next hour of the day and the the errors to get distribution
univariate_forecast_t[i + 1, ] <- predict_point[i + 1] + mu_sigma[i + 1, "mu"] + inverse_ecdf(quantiles, i) * mu_sigma[i + 1, "sigma"]
}
# Dependence Learning: Pairing up the Copula
multivariate_forecast <- matrix(nrow = nrow(R_emp), ncol = ncol(R_emp))
for (i in seq(1:24)) {
j <- 1
for (r in R_emp[, i]) {
multivariate_forecast[j, i] <- univariate_forecast_t[i, r]
j <- j + 1
}
}
print("Prediction DONE...")
# Extracting the peaks from the multivariate distribution
peaks_dis_B <- apply(multivariate_forecast, MARGIN = 1, FUN = max)
# Extracting the peak from the test day
peak <- max(df_test[, "y"])
hist(peaks_dis_B, main = "Histogram of Peaks Distribution B", xlab = "Peaks", breaks = "Sturges")
return(crps_sample(peak, peaks_dis_B))
}
# specify the length for rolling iterations in days
len_test <- 28
crps_scores <- list()
for (d in seq(1, len_test)) {
crps_score <- getCRPS_B(d)
# Append the CRPS score to the list
crps_scores[[d]] <- crps_score
}
## [1] 1
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 2
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 3
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 4
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 5
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 6
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 7
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 8
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 9
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 10
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 11
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 12
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 13
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 14
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 15
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 16
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 17
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 18
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 19
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 20
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 21
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 22
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 23
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 24
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 25
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 26
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 27
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 28
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
print("Mean CRPS for 28 days in 2024 for NeuralProphet")
## [1] "Mean CRPS for 28 days in 2024 for NeuralProphet"
print(mean(unlist(crps_scores)))
## [1] 2916.514
library(tidyverse)
library(scoringRules)
library(ggplot2)
library(copula)
# This is a first general implementation of version B
# Load the model data
# Here: NeuralProphet
# model <- read.csv("../data/forecasts/load_22-24_model-neuralprophet_2024IsForecasted.csv")
# model$ds <- as.POSIXct(model$ds, tz = "UTC")
# Here: Arma
model <- read.csv("../data/forecasts/load_22-24_model-arma.csv")
model$ds <- as.POSIXct(model$ds, tz = "UTC")
# specify the length for the error learning phase in days
length <- 365
# Filter only for the years 2023 (Training) and 2024 (Testing)
model <- model |>
filter((year(ds) == 2023) | (year(ds) == 2024))
getCRPS_B <- function(d) {
print(d)
# --------- Error Learning Phase ---------
# Getting the test data: starting from day i get the next 365 days
# Achtung: Hier dürfen nur Training oder Testing data verwendet werden
df_train <- model[((d - 1) * 24 + 1):(((364 + d) * 24) - 1), ]
mu_sigma <- data.frame(hour = seq(0, 23), mu = rep(0, 24), sigma = rep(0, 24))
for (i in seq(0, 23)) {
load_h <- df_train |>
filter(hour(ds) == i)
mean <- mean(load_h$residuals)
std <- sqrt(var(load_h$residuals))
mu_sigma[i + 1, "mu"] <- mean
mu_sigma[i + 1, "sigma"] <- std
# Extract the residuals from the data
resid <- data.frame(resid = load_h$residuals)
# Standardize the residuals
resid <- resid |>
mutate(resid_std = (resid - mu_sigma[i + 1, "mu"]) / mu_sigma[i + 1, "sigma"])
# Generate the ECDF
assign(paste0("ecdf_", i), ecdf(resid$resid_std))
}
# Generate the Inverse ECDF
# Input: p... Probabilty , i... hour
# Output: quantile
inverse_ecdf <- function(p, i) {
quantile(get(paste0("ecdf_", i)), p, names = FALSE)
}
print("Error Learning DONE...")
# ---- Dependence Learning Phase ----
# Define the length m of the learning phase
m <- 90
# Initialize the training sample of the past m errors for every hour to learn the copula
X <- matrix(nrow = m, ncol = 0)
for (i in seq(0, 23)) {
# get the hour data
load_h <- df_train |>
filter(hour(ds) == i)
# get the errors
epsilon <- as.vector(load_h$residuals)
# Select only the last m errors
epsilon <- epsilon[(length(epsilon) - m):length(epsilon)]
# standardize the errors
epsilon_std <- (epsilon - mu_sigma[i + 1, "mu"]) / mu_sigma[i + 1, "sigma"]
# Apply the ecdf
ecdf <- get(paste0("ecdf_", i))
epsilon <- ecdf(epsilon_std)
X <- cbind(X, epsilon)
}
# Create the Copula and the Rank Matrix
# 2) Non-Parametric Approach
# Consider X as the sample from the empirical copula
# Apply the rank function to get the rank matrix
R_emp <- apply(X, 2, rank)
print("Dependence Learning DONE...")
# ------- Prediction phase -------
# Generate the 24 univariate for every hour
# m also defines the quantile level i.. 1 - m with i/m+1 quantiles
univariate_forecast_t <- matrix(nrow = 24, ncol = m)
quantiles <- seq(1 / (m + 1), m / (m + 1), 1 / (m + 1))
# Get the next 24 hours after the testing period
# Hier nur Testing data vom Model!!!
df_test <- model[((d + length - 1) * 24):((d + length) * 24 - 1), ]
predict_point <- df_test[, "yhat"]
for (i in seq(0, 23)) {
# Get the prediction for the next hour of the day and the the errors to get distribution
univariate_forecast_t[i + 1, ] <- predict_point[i + 1] + mu_sigma[i + 1, "mu"] + inverse_ecdf(quantiles, i) * mu_sigma[i + 1, "sigma"]
}
# Dependence Learning: Pairing up the Copula
multivariate_forecast <- matrix(nrow = nrow(R_emp), ncol = ncol(R_emp))
for (i in seq(1:24)) {
j <- 1
for (r in R_emp[, i]) {
multivariate_forecast[j, i] <- univariate_forecast_t[i, r]
j <- j + 1
}
}
print("Prediction DONE...")
# Extracting the peaks from the multivariate distribution
peaks_dis_B <- apply(multivariate_forecast, MARGIN = 1, FUN = max)
# Extracting the peak from the test day
peak <- max(df_test[, "y"])
hist(peaks_dis_B, main = "Histogram of Peaks Distribution B", xlab = "Peaks", breaks = "Sturges")
return(crps_sample(peak, peaks_dis_B))
}
# specify the length for rolling iterations in days
len_test <- 28
crps_scores <- list()
for (d in seq(1, len_test)) {
crps_score <- getCRPS_B(d)
# Append the CRPS score to the list
crps_scores[[d]] <- crps_score
}
## [1] 1
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 2
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 3
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 4
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 5
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 6
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 7
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 8
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 9
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 10
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 11
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 12
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 13
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 14
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 15
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 16
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 17
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 18
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 19
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 20
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 21
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 22
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 23
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 24
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 25
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 26
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 27
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
## [1] 28
## [1] "Error Learning DONE..."
## [1] "Dependence Learning DONE..."
## [1] "Prediction DONE..."
print("Mean CRPS for 28 days in 2024 for ARMA")
## [1] "Mean CRPS for 28 days in 2024 for ARMA"
print(mean(unlist(crps_scores)))
## [1] 16453.73
